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Abstract. A rigorous homogenization theory of metamaterials - artificial periodic structures 
judiciously designed to control the propagation of electromagnetic waves - is developed. All coarse- 
grained fields are unambiguously defined and effective parameters are then derived without any 
heuristic assumptions. The theory is an amalgamation of two concepts: Smith & Pendry's physical 
insight into field averaging and the mathematical framework of Whitney-Nedelec-Bossavit-Kotiuga 
interpolation. All coarse-grained fields are defined via Whitney forms and satisfy Maxwell's equa- 
tions exactly. The new approach is illustrated with several analytical and numerical examples and 
agrees well with the established results (e.g. the Maxwell-Garnett formula and the zero cell-size 
limit) within the range of applicability of the latter. The sources of approximation error and the 
respective suitable error indicators are clearly identified, along with systematic routes for improving 
the accuracy further. The proposed approach should be applicable in areas beyond metamaterials 
and electromagnetic waves - e.g. in acoustics and elasticity. 

1 Introduction 

Electromagnetic and optical metamaterials are periodic structures with features smaller than the 
vacuum wavelength, judiciously designed to control the propagation of waves. Typically, resonance 
elements (variations of split-ring resonators) are included to produce nontrivial and intriguing effects 
such as backward waves and negative refraction, cloaking, slow light ( "electromagnetically induced 
transparency"), and more; see [2, 8, 11, 20, 23, 29, 38, 39, 40], to name just a few representative 
publications, and references there. 

To gain insight into the behavior of such artificial structures and to be able to design useful 
devices, one needs to approximate a given metamaterial by an effective medium with dielectric 
permittivity e c g and magnetic permeability // e g- or, in more general cases where magnetoelectric 
coupling may exist, with a 6 x 6 parameter matrix. A variety of approaches have been explored 
in the literature (e.g. [22, 24, 25, 26, 28, 31], again to name just a few publications), with notable 
accomplishments in designing novel and interesting devices and structures, as cited above. However, 
the existing approaches are mostly heuristic, and there is a clear need for a consistent and rigorous 
theory - rigorous in the sense that all "macroscopic" (coarse-grained) fields are unambiguously and 
precisely defined, giving rise to equally well defined effective parameters. 
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It is the main objective of this paper to put forward such a theory. The methodology advocated 
here is an amalgamation of two very different lines of thinking: one relatively new and driven 
primarily by physical insight, and the other one well established and mathematically rigorous. The 
physical insight is due to Pendry & Smith [28], who prescribed different averaging procedures for 
the h and b fields (and similarly for the e and d fields). The practical results have been excellent, 
and yet it has remained a bit of a mystery why, say, h and b must be averaged differently. The 
justification for that in [28] and other publications comes from the analogy with staggered grid 
approximations in finite difference methods, but it is unclear why physics should be subordinate to 
numerical methods and not the other way around. 

Remark. Throughout the paper, small letters b,e,d, etc., will denote the "microscopic" - i.e. 
true physical - fields that in general vary rapidly as a function of coordinates. Capital letters will 
refer to smoother fields, to be defined precisely later, that vary on the scale coarser than the lattice 
cell size. 

The second root of the proposed methodology is the mathematical framework developed by 
Whitney in the middle of the 20th century [37]. A few decades later Whitney's theory was discovered 
to be highly relevant in computational electromagnetism, due primarily to the work of Nedelec, 
Bossavit and Kotiuga [6, 7, 12, 18, 19]. It should be emphasized, however, that in the present 
paper the Whitney-Nedelec-Bossavit-Kotiuga (WNBK) framework is used not for computational 
purposes but to define, analytically, the coarse-grained fields. 

The end result of combining WNBK interpolation with Smith & Pendry's insight is a mathemat- 
ically and physically consistent model that is rigorous, general (e.g., applicable to magnetoelectric 
coupling) and yet simple enough to be practical. The theory is supported by a number of an- 
alytical and numerical case studies and is consistent with the existing theories and results (e.g. 
with the Maxwell-Garnett mixing formula and with the zero cell-size limit) within the ranges of 
applicability of the latter. Approximations that have to be made, and the respective sources of 
error, are clearly identified. Several routes for further accuracy improvement are apparent from the 
theoretical analysis. 

2 Some Pitfalls 

Prior to discussing what needs to be done, it is instructive to review what not to do. Some 
averaging procedures appear to be quite natural and yet upon closer inspection turn out to be 
flawed. Physically valid alternatives are introduced in the subsequent sections. 

Passing to the limit of the zero cell size. A distinguishing feature of the homogenization 
problem for metamaterials is that the cell size a, while smaller than the vacuum wavelength Ao at a 
given frequency u, is not vanishingly small. A typical range in practice is a ~ 0.1 — 0.3Ao- Therefore 
classical homogenization procedures valid for a — > - e.g. Fourier [27] or two-scale analysis [3] - 
have limited or no applicability here. Independent physical [16] and mathematical [34] arguments 
show that the finite cell size is a principal limitation rather than just a constraint of fabrication. 
If the cell size is reduced relative to the vacuum wavelength, the nontrivial physical effects (e.g. 
"artificial magnetism") ultimately disappear, provided that the intrinsic dielectric permittivity e of 
the materials remains unchanged. On physical grounds [34], this can be explained by the operating 
point on the normalized Bloch band diagram falling on the "uninteresting" acoustic branch. 

Mollifiers. The classical approach to defining the macroscopic fields is via convolution with a 
smooth mollifier function (e.g. Gaussian-like) [21]. This necessitates an intermediate scale for the 
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Figure 1: The field and magnetization in the direction normal to the interface. 

mollifier, much coarser than the cell size but still much finer than the wavelength in the material. 
For natural materials, this requirement is easily fulfilled because the cell size is on the order of 
molecular dimensions and much smaller than the wavelength. In contrast, for metamaterials the 
cell size is typically ~0.1 - 0.3 of the wavelength, and no intermediate scale is available for the 
mollifier. 

Averaging over the cell. The following argument shows that simple cell- averaging of the 
fields is for metamaterials inadequate. To fix ideas, consider a simplified picture first. A qualitative 
behavior of a tangential component of the microscopic b-field in the direction normal to the surface 
is shown in Fig. 1. For our purposes at this point, the precise field distribution is unimportant; one 
may have in mind, for example, a single Bloch wave moving away from the surface, even though 
later on it will be important to consider a superposition of Bloch waves. 

Clearly, the average field Bo = (b) over the cell is not in general equal to the field b(0+) inside 
the material immediately at the boundary. Yet it is b(0+) that couples with the field in the air: 
b(0+) = b(0— ), where b(0— ) is the field in the air immediately at the boundary. (Intrinsically 
nonmagnetic materials are assumed throughout.) The classical boundary condition is recovered 
if an auxiliary H field is introduced in such a way that H(0+) = b(0+); magnetization then is 
47rm = b — H, as schematically shown in the figure. In other words, it is the difference between 
the cell-averaged field b and its value at the boundary that ultimately defines magnetization and 
the permeability. 

This picture, albeit simplified (in particular, it ignores a complicated surface wave at the inter- 
face [25, 14]), does serve as a useful starting point for a proper physical definition of field averages 
and effective material parameters. 

Since the cell-averaged field is, for a non-vanishing cell size, not generally equal to its boundary 
value and does not satisfy the proper continuity condition at the interface, the use of such fields in a 
model would result in nonphysical equivalent electric / magnetic currents on the surface (jumps of a 
tangential component of the magnetic / electric field), nonphysical electric/magnetic surface charges 
(jumps of the normal component of D or B), and incorrect reflection/transmission conditions at the 
boundary. It is true that, as a "zero-order" approximation, the cell-averaged field is approximately 
equal to its pointwise value; however, equating them means ignoring the variation of the fields over 
the cell and hence throwing away the very physical effects under investigation. 



"Magnetic dipole moment per unit volume." This textbook definition of magnetization 
turns out to be flawed as well. Simovski & Tretyakov [26] give a counterexample for a system 
of two small particles, but in fact their argument is general. Suppose that a large volume of a 
metamaterial has been in some way homogenized and is now represented, to an acceptable level of 
approximation, by effective parameters /i e g, e e ff. Consider then a standing electromagnetic wave 
in this material (as produced e.g. in a cavity or by reflection off a mirror). It is well known 
that the electric and magnetic fields in a lossless standing wave are shifted by one quarter of the 
wavelength. At a node of the electric field, i.e. at a point in space where the electric field is zero, 
the magnetic flux density and hence the magnetization are maximum. Yet the zero electric field at 
the node implies zero currents (both polarization and conduction) and hence zero magnetic dipole 
moment (as there is no spin-related intrinsic moment by assumption). This inconsistency shows 
that magnetization cannot be defined as the dipole moment per unit volume. 

Bulk parameters. It is known that even for a homogeneous isotropic infinite medium the pair 
of parameters e and \i are not defined uniquely. Indeed, the total microscopic current j can be split 
up - in principle, fairly arbitrarily - into the "electric" and "magnetic" parts, j = dtp + cV x m 
[36, 15, 14]. The h field is then defined accordingly, as b — 47rm, giving rise to the respective value of 
/i e fj that depends on the choice of m. A more general "Serdyukov-Fedorov" transformation leaves 
Maxwell's equations invariant but changes the values of the material parameters [4, 5, 36]: 

d' = d + VxQ, h' = h + c _1 a t Q 

b' = b + V x F, e' = e + cT^F 

where Q and F are arbitrary fields (with a valid curl). It is possible [1, 36] to set [i = 1, in which 
case the dielectric function carries all relevant information but becomes spatially dispersive (its 
transform depends on the spatial Fourier vector k). 

Thus even for a homogeneous infinite medium it is only the product e/i that is unambigu- 
ously defined, with its direct physical relation to phase velocity v p = */ejl. The situation changes 
thoroughly when a material interface (for simplicity, with air) is considered. Classical boundary 
conditions for the tangential continuity of the H and E fields 1 fix the ratio of the material param- 
eters via the intrinsic impedance rj = \J Jt/e, which, taken together with their product, identifies 
these two parameters separately and uniquely. 

It is clear, then, that any complete and rigorous definition of the effective electromagnetic 
parameters of metamaterials must account for boundary effects. 

3 Coarse- Grained Fields 
3.1 The Guiding Principles 

Consider a periodic structure composed of materials that are assumed to be (i) intrinsically non- 
magnetic (which is true at sufficiently high frequencies [13, 16]); (ii) satisfy a linear local constitutive 
relation d = ee. For simplicity, we assume a cubic lattice with cells of size a. 

Maxwell's equations for the microscopic fields are, in the frequency domain and with the 
exp(— iuit) phasor convention, 

V x e = \ujc~ l h 

1 These conditions are local and much simpler than the ones arising in the model with fj, = 1 [1, 36]. 
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V x b = — \ujc d 

The coarse-grained fields B, H, E, D must be defined in such a way that the boundary conditions 
are honored. (As already noted, simple cell-averaging does not satisfy this condition.) From the 
mathematical perspective, these fields must lie in their respective functional spaces 

E,He iJ(curl, fl); B, D G H(div, Q) (1) 

where O is a domain of interest that for mathematical simplicity is assumed finite. Symbol will 
henceforth be dropped to shorten the notation. Rigorous definitions of these functional spaces are 
available in the mathematical literature (e.g. [17]). From the physical perspective, constraints (1) 
mean that the E, H fields possess a valid curl as a regular function (not as a Schwartz distribution), 
while B and D have a valid divergence. This implies, most importantly, tangential continuity of 
E, H and normal continuity of B and D across material interfaces. The fields in i?(curl) are also 
said to be curl- conforming, and those in H(div) to be div- conforming. 

In principle, any choice of a curl-conforming H field produces the respective "magnetization" 
47rm = b — H and leaves Maxwell's equations intact. However, most of such choices will result 
in technically valid but completely impractical and arbitrary constitutive laws, with the "material 
parameters" depending more on the choice of H than on the material itself. 

Below, I argue that construction of the coarse-grained fields via the Whitney-Nedelec-Bossavit- 
Kotiuga (WNBK) interpolation has particular mathematical and physical elegance, which leads to 
practical advantages in the computation of fields in periodic structures. 

3.2 Background: WNBK Interpolation 

For a rigorous definition of the coarse-grained H-field, we shall need an interpolatory structure 
referred to in the literature as the "Whitney complex" [6, 7]; however, acknowledgment of the 
seminal contributions of Nedelec, Bossavit and Kotiuga [18, 6, 12] is quite appropriate and long 
overdue. Still, the subscripts of quantities related to the Whitney complex will for brevity be just 
'W rather than 'WNBK'. 

WNBK complexes form a basis of modern finite element analysis with edge and facet elements. 
However, our objective here is not to develop a numerical procedure; rather, the mathematical 
structure that has served so well in numerical analysis is borrowed and applied to fields in meta- 
material cells. 

The original Whitney forms [37] are rooted very deeply in differential geometry and algebraic 
topology, and the interested reader can find a complete mathematical exposition in the literature 
cited above. For the purposes of the present paper, we need a small subset of this theory where 
the usual framework of vector calculus is sufficient. It will also suffice to consider a reference cell 
as a cube [— 1,+1] 3 . This shape can be transformed into an arbitrary rectangular parallelepiped 
by simple scaling and, if necessary, to any hexahedron by a linear transformation as in [35]. 

We shall need two interpolation procedures: (i) given the circulations [q] a = f a qdl of a field q 
over the 12 edges of the cell (a = 1, 2, ... , 12), extend that field into the volume of the cell; and (ii) 
given the fluxes [[q]]p = fp q • dS of a field q over the six faces of the cell [fi = 1, 2, . . . , 6), extend 
that field into the volume of the cell. Single brackets denote line integrals over an edge; double 
brackets - surface integrals over the faces. Typically, item (i) will apply to the e and h fields and 
(ii) - to the d and b fields. 
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Figure 2: A 2D analog of the vectorial interpolation function w Q (in this case, associated with the 
central vertical edge shared by two adjacent cells). Tangential continuity of this function is evident 
from the arrow plot; its circulation is equal to one over the central edge and to zero over all other 
edges. 



Consider an edge a along a ^-direction, where £ is one of the coordinates x, y, z, and let r\ and 
r be the other two coordinates, with (£,77, r) being a cyclic permutation of (x,y,z) and hence a 
right-handed system. The edge a is then formally defined as 

-1<£<1; r] = r] a ; r = T a ; r] a ,T a = ±l 

Associated with this edge is a vectorial interpolating function 

w a = - (l + r/ a r/)(l + r Q r)e (2) 

where the hat symbol denotes the unit vector in a given direction. For illustration, a 2D analog of 
this vector function is shown in Fig. 2. 

In 3D, there are 12 of such interpolating functions - one per edge - in the cell. It is straight- 
forward to verify that the edge circulations of these functions have the Kronecker-delta property 

[w a ] a > = 5 aa < (3) 

Indeed, on edge a itself, r] a = r a = r\ = r = 1; hence the value of w a on its "own" edge is which 
produces unit circulation upon integration over the edge of length two ([— 1,+1]). For edges not 
along the £ direction, the circulation is obviously zero because the vector function w Q is orthogonal 
to them. For other edges in the £ direction, one of the coordinates rj or r is equal to —1 and the 
respective factor in (3) evaluates to zero, again producing a zero circulation. 
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Figure 3: A 2D analog of the vectorial interpolation function (in this case, associated with the 
central vertical edge). Normal continuity of this function is evident from the arrow plot; its flux is 
equal to one over the central edge and zero over all other edges. 

The Kronecker-delta property guarantees that the interpolating functions are linearly inde- 
pendent over the cell and span a 12-dimensional space of vectors that can all be represented by 
interpolation from the edges into the volume of the cell: 



We shall call this 12-dimensional space W cur i: the 'W honors Whitney and 'curl' indicates fields 
whose curl is a regular function rather than a general distribution. This implies, in physical terms, 
the absence of equivalent surface currents and the tangential continuity of the fields involved. Any 
adjacent lattice cells sharing a common edge will also share, by construction of interpolation (4), 
the field circulation over that edge, which ensures the continuity of the tangential component of 
the field across all edges. 

The curls of w a are not linearly independent but rather, as can be demonstrated, lie in the 
six-dimensional space spanned by the following functions: 
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A 2D analog of a typical function v is shown in Fig. 3. 

The v functions satisfy the Kronecker delta property with respect to the face fluxes: 




(6) 



Indeed, the flux of vi = x(l + x)/2 is obviously equal to one over the face x = 1 and zero over 
all other faces; similar conditions hold for the other five functions. Hence one can consider vector 



interpolation from fluxes on the cell faces into the volume of the cell, conceptually quite similar to 
edge interpolation (4): 

6 

q = E [WW (7) 

The six-dimensional space spanned in a lattice cell by the v functions will be denoted with Wdi v > 
reflecting the easily verifiable fact that the normal component of the field interpolated via (7) is 
continuous across the common face of two adjacent cells; hence the generalized divergence of this 
field exists as a regular function, not just as a distribution (that is, there are no (5-sources - surface 
charges - on the cell faces). 

Importantly, as already mentioned, the curls of functions from W CUT \ lie in Wdi v , or symbolically 
in terms of the functional spaces, 

V x W curl E W div (8) 
To summarize, the following properties are critical for our construction of the coarse-grained fields: 

1. Twelve functions w (one per cell edge) interpolate the field from its circulations over the 
edges into the cell. The resulting field is tangentially continuous across all cell boundaries. 
The w functions span a 12-dimensional functional space W cur i. 

2. Six functions v (one per cell face) interpolate the field from its fluxes through the faces into 
the cell. The resulting field has normal continuity across all cell faces. The v functions span 
a 6-dimensional functional space Wdiv 

3. Any vector field in W cur i is uniquely defined by its twelve edge circulations. 

4. Any vector field in Wdw is uniquely defined by its six face fluxes. 

5. V x W cur i E Wd iv . 

3.3 Construction of the Coarse- Grained Fields 

As previously noted, naive averaging of the E and H fields over the cell adjacent to a material- 
air interface breaks the tangential continuity of these fields across the boundary and is therefore 
unacceptable. WNBK interpolation provides a suitable alternative. Let us define the coarse-grained 
fields as WNBK interpolants of the actual edge circulations via the w functions, in accordance with 
(4). Within each lattice cell, 

12 

E = £[e] a w a = Wcurl([e]l-12) (9) 
a=l 

with a completely similar expression for the H- field. The '=' signs indicate that this is the definition 
of E and H, as well as of the WNBK curl- interpolation operator W CU ri- 

Similarly, the B and D fields are defined as interpolants in Wdi v ; this interpolation, from the 
actual face fluxes into the cells, is effected by the v functions; within each cell, 

6 

B = EtMW = WdivdMh-e) (10) 
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and a completely analogous expression for the D field. 

We now decompose the microscopic fields e, b, d into coarse-grained parts E, B, D defined 
above and rapidly varying remainders e~, b~, d~: 

e = E + e~ (11) 

d = D + d~ (12) 
b = B + b~ (13) 
Importantly, b has an alternative decomposition where H is taken as a basis: 

b = H + 4vrm~ (14) 

With these splittings, Maxwell's equations become 

Vx(E + e~) = iw(B + b~) (15) 

Vx(H + 47rm~) = - iu/(D + d~) (16) 

It is at this point that the role of the WNBK interpolation becomes apparent: the scales separate. 
Indeed, E is, by construction, in W cm i, and therefore V x E is in Wdi v , and so is, by construction, 
B. In that sense, the capital-letter terms in (15) are fully compatible. Furthermore, E - again by 
definition - has the same edge circulations as the original microscopic field e; hence, for any face 
of any cell, 

f (V x E) • dS = f E • dl = \ujc~ 1 J b • dS 

J face J face edges J face 

where the Stokes theorem and the microscopic Maxwell equation were used. But the B-field has the 
same face fluxes as b by construction and, since these face fluxes define the field in Wdiv uniquely, 
we have 

V x E = io;c _1 B (17) 

Analogously, 

V x H = - iwc _1 D (18) 

Thus the coarse level has separated out and, remarkably, the WNBK fields satisfy the Maxwell 
equations, as well as the proper continuity conditions, exactly. The underlying reason for that is 
the compatibility of the curl- and div-interpolations, i.e. condition (8). 

For the rapidly changing components, straightforward algebra yields, from (15) and (16): 

V x e~ = kdb~ (19) 
4fVxm~ = -iwd~ (20) 

with the "constitutive relationships" 

d~ = ee~ + (eE - D) (21) 

4vr m~ = b~ + (B - H) (22) 
All edge circulations of e~, m~ and all face fluxes of b~, d~ are zero. For example, 

[e~] = for each edge (23) 
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Figure 4: The lattice (with inclusions of arbitrary shape) serves as a "scaffolding" for the construc- 
tion of coarse-grained fields. The curl-conforming fields such as H are interpolated into the cells 
from the edge circulations, while the div-conforming fields such as B are interpolated from the face 
fluxes. 

If/when the coarse-grained fields have been found from their Maxwell equations, one may convert 
the two equations for the rapid fields into a single equation for e~: 

V x V x e~ - w 2 ee~ = w 2 (eE - D) - iwV x (B - H) (24) 

This equation can be solved for each cell, with the Dirichlet-type zero-circulation boundary condi- 
tions (22). In principle, this fast-field correction will increase the accuracy of the overall solution. 
However, the remainder of the paper will be focused on our main subject - the coarse fields and 
the corresponding effective parameters. 

4 Material Parameters 

4.1 Procedure for the Constitutive Matrix 

We are now in a position to define effective material parameters, i.e. linear relationships between 
D,B and E,H. Consider a metamaterial lattice with inclusions of arbitrary shape. We construct 
a coarse-grained curl-conforming field such as H using the lattice cell edges as a "scaffolding": 
the field is interpolated into the cells from the edge circulations of the respective microscopic field 
(Fig. 4). Thus, the edge circulations of the microscopic and the coarse-grained curl-conforming 
fields are the same. Similar considerations apply to div-conforming fields, but with interpolation 
from the faces. 

Let the electromagnetic field be approximated, in a certain region, as a linear combination of 
some basis waves ip a : 

a a 
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In the most general case, and all ip a are six-component vector comprising both microscpoic 
fields; e.g. ^ eh = {ip e ,ip h }, etc. However, in the absence of magnetoelectric coupling it is natural 
to consider each pair of fields (e, d) and (h, b) separately, and then the ips have three components 
rather than six. The basis waves ip a can, but do not have to, be Bloch waves. While Bloch waves 
are most general, using multipole expansions with respect to a central inclusion in the lattice cell 
could be beneficial in some cases. 

To each basis wave, there corresponds a WNBK interpolant described above. In the simple case 
of the H and B fields aligned in the same direction, the pointiwse material parameters are found 
as the ratios of B and H and of D and E. 

To define the material parameters in a more general setting, the field ratios above need to 
be replaced with "generalized ratios," as follows. At any given point r in space, consider the 
WNBK curl-interpolants E a (r) = W CUT \([ip%]i-i2) and H a (r) = W cur i([V>a]i-i2) for each basis 
wave ip a . Similarly, consider the WNBK div-interpolants D a (r) = Wdi v ([[da]]i-6) and B a (r) = 
Wdi V ([[b a ]]i_ fl ). 

Then, for each basis wave a, we seek a linear relation 

^ B (r) = C(r)V^(r) 

where £ is a 6 x 6 constitutive matrix characterizing, in general, anisotropic material behavior with 
(if the off-diagonal blocks are nonzero) magnetoelectric coupling. Similar to tp db , the six-component 
vectors ip DB comprise both fields, but coarse-grained; same for ijj EH . In matrix form, the above 
equations are 

^ DB (r) = C(r)^ EH {r) (25) 

where each column of the matrices ^ DB and ^> EH contains the respective basis function. (Illustra- 
tive examples in subsequent sections may help to clarify these notions and notation.) 

If exactly six basis functions are chosen, one obtains the constitutive matrix by straightforward 
matrix inversion; if the number of functions is more than six, the pseudoinverse [10] is appropriate: 

£(r) = y DB (r)(V EH ) + (r) (26) 

The approximate parameters for each cell are then found simply by cell- averaging of C( r )- Quite a 
notable feature of this construction is that the coarse-grained fields corresponding to a particular 
basis wave satisfy the Maxwell equations with this material parameter £(r) exactly. 

One may wonder why such cell-averaging could not be done in the very beginning, for the 
original "microscopic" Maxwell equations - in which case one would get the trivial value of unity 
by averaging the intrinsic permeability ji = 1. To see the difference, assume for simplicity that 
there is no magnetoelectric coupling (this assumption does not change the essence of the argument) 
and compare the curl-curl equations for the microscopic and coarse-grained fields: 



V x V x e - 




V x C^'WV x E - (~) C e (r)E = 

where Q e and are the electric and magnetic submatrices of the £(r) defined above. The micro- 
scopic equation can indeed be averaged over the cell (noting that the curl and averaging operators 
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commute), but this would lead to the cell-averaged field that for metamaterials is inadequate, as 
has already been discussed. 

In the case of the coarse-grained field, the splitting of the material matrix into its cell average 
plus a fluctuating component, C = (C) + C~i does make sense. A perturbation analysis shows that 
this splitting produces fictitious sources like £~(r)E that approximately average to zero within each 
cell because the coarse-grained field varies slowly and the mean of £~ is zero by definition. The 
field due to these spurious equivalent sources is therefore small. This qualitative argument is not a 
complete substitute for a detailed mathematical analysis that could be undertaken in the future. 

4.2 The Recipe 

Let the vectorial dimension of the problem - i.e. the total number of vector field components - be 
N. Most generally in electromagnetism, N = 6 (three components of E and three of H), but if 
only one field is involved, then N = 3, and if that field has only one component, then N = 1, etc. 
Let M be the total number of basis waves; since M < N makes no practical sense, we shall always 
assume M > N. 

It is convenient to summarize the procedure described above as a "recipe" for finding the 
effective parameters: 

• Choose a set of basis waves ip a that provide a good approximation of the fields within a cell 
well. In general, Bloch waves are suitable candidates, although other choices may be possible 
under specific circumstances. The number of basis functions should be equal to or greater 
than the vectorial dimension N of the problem. 

• Find the curl- and div-WNBK interpolants of each basis wave. (This step requires the com- 
putation of face fluxes and edge circulations of the respective fields in the wave.) 

• Assemble these WNBK interpolants into the ^ DB and ty EH matrices of (25). 

• Find the coordinate-dependent parameter matrix from (26). 

• The cell-average of this matrix gives the final result: an N x N (6 x 6 in the most general 
case) matrix of effective parameters. 

4.3 Errors and Error Indicators 

Following the analysis of the previous sections, one can identify the sources of error in the re- 
placement of the actual metamaterial with an effective medium. Even more importantly, ways to 
improve the accuracy can also be identified, as well as the limitations of such improvement. In 
the remainder, we shall ignore the usual numerical errors of evaluating the basis waves (e.g. finite 
element or finite difference discretization errors if these methods are used to compute the basis) be- 
cause such errors are already understood quite well and are only tangentially related to the subject 
of this paper. 

Three sources of error can be distinguished in the proposed homogenization procedure. 

1. In-the-basis error. If the number of basis waves is strictly greater than the vectorial dimension 
of the problem (M > N), then system (25) does not generally have an exact solution and is 
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solved in the least-squares sense, as the pseudoinverse in (26) indicates. From (25) and (26), 
the discrepancy between the fields is (with the dependence on r implied) 

where / is the 6x6 identity matrix. Therefore the matrix norm \\I — (\& EH }+\& EH \\ is a 
suitable indicator of the in-the-basis error. If M = N (for example, six basis waves in a 
generic problem of electrodynamics), then this matrix is normally zero and the in-the-basis 
error vanishes. 

2. Out-of-the-basis error. Any field can be represented as a linear combination of the basis 
waves, plus a residual field. If a good basis set is chosen, this residual term, and hence the 
error associated with it, is small. Any expansion of the basis carries a trade-off: the residual 
field and the out-of-the-basis error will decrease, but the in-the-basis error may increase. 

3. Parameter averaging. As an intermediate step, the proposed homogenization procedure yields 
a coordinate-dependent parameter matrix £(r) that is in the end averaged over the lattice 
cell. This is discussed in the end of the previous section. 

The limitations of the effective-medium approximation now become apparent as well. If a 
sufficient number of modes with substantially different characteristics exist in a metamaterial (e.g. 
in cases of strong anisotropy), then the in-the-basis error will be high. This is not a limitation of 
the specific procedure advocated in the paper, but rather a reflection of the fact that the behavior 
of fields in the material is in such a case too rich to be adequately described by a single effective- 
parameter tensor. 

On a positive side, several specific ways of improving the approximation accuracy can be iden- 
tified. Obviously, no accuracy gain is completely free; better approximations may require more 
information than a single material matrix can contain. 

• The cell problem (19)- (22) for the rapidly varying fields can be solved, once the coarse-grained 
fields have been found from the effective-parameter model. 

• If the relative weights of different modes in a particular model can be estimated a priori at 
least roughly, then system (25) can be biased toward these modes. The downside of such 
biasing is that the material parameters are no longer a property of the material alone but 
become partly problem-dependent. 

• The size and composition of the basis set can be optimized for common classes of problems, 
once practical experience of applying the methodology of this paper is accumulated. 

• The last step of the proposed procedure - the cell-averaging of matrix £ - can be omitted. 
The problem for the coarse-grained fields is then exact. The trade-off is that the numerical 
solution of such a problem may be computationally expensive, as the material parameter 
varies within the cell. There is of course a spectrum of practical compromises where £ would 
be approximated within a cell not to order zero (i.e. as a constant) but to some higher order. 

• One may envision adaptive procedures whereby the basis waves are updated after the problem 
has been solved, and then new material parameters are derived from the new basis set. 
Recursive application of adaptivity may result in a substantial improvement of the solution. 
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The downside, in addition to the increased cost, is the same as before: the material parameters 
become partly problem-dependent. 



In connection with the last item, it is clear that the basis set should in general reflect the 
symmetry and reciprocity properties of the problem. In particular, the basis should as a rule 
include pairs of waves traveling in the opposite directions. 

4.4 Causality and Passivity 

Physical considerations indicate that the proposed procedure should be expected to produce causal 
and passive effective media, at least if M = N. Indeed, suppose that the opposite is true - e.g. 
£g ff < (with anisotropy neglected for simplicity), violating passivity. The effective parameters, 
however, apply by construction to all basis waves exactly; hence e" s < would imply power 
generation in the actual physical modes in the actual passive metamaterial, which is impossible. 

While such physical considerations are quite plausible, a rigorous mathematical analysis is still 
desirable in future research. 

5 Verification 

5.1 Empty Cell 

This first test can be viewed as a "sanity check" : will the proposed procedure produce unit perme- 
ability and permittivity for an empty cell? This is not a completely trivial question: it is common 
for the alternative methods in the literature to yield spurious Bloch-like factors that are then 
removed from the effective parameters by fiat. 

For a generic plane wave passing through an empty cell in either 2D or 3D, direct calculation 
that follows the proposed methodology shows that exact parameters /x e g = 1 e c fj = 1, without any 
spurious factors and of course with no magnetoelectric coupling, are indeed obtained, regardless 
of the direction of the plane wave. No fudge factors or heuristic adjustments are needed to bring 
these parameters to unity. 

5.2 Example: One-component static fields 

This is an obvious case that serves just as an illustrative example and a consistency check for the 
proposed methodology. Let a static field (for instance, electrostatic) have only one component 
(say, z) that must be independent of z due to the zero-divergence condition. Then the lattice 
cells become effectively two-dimensional and the div-conforming WNBK interpolant for d reduces 
just to a constant Do whose flux through the cell is equal to the flux of d. The curl-conforming 
WNBK interpolant W COT ./(E) would generally be a bilinear function of x and y, but since the field 
is constant in the xy-plane, this interpolant also reduces to a constant Eo = E. 
The dielectric permittivity thus is 




SE 




exactly as should be expected. 
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5.3 Example: Comparison with the Maxwell-Garnett formula 

The next test is consistency with the classical Maxwell-Garnett (M-G) mixing formula for a two- 
component medium, in the limit where this formula is valid. Namely, the fill factor for the "inclu- 
sion" component in a host material is assumed to be small and the field is static. (We shall not deal 
with radiative corrections to the polarizability and to M-G in this paper.) The M-G expression for 
the effective permittivity is, in 2D, 

i + fx 
1 - fx 

where x 1S the polarizability of inclusions in a host medium. In particular, for cylindrical inclusions 
in a non-polarizable host, x = ( e cyi — 1)/ ( e cyi + 1) an d the M-G formula becomes 

ecyi(l + /) + (!-/) 
CMG ' 2D e cyl (l -/) + (! + /) 

The proposed methodology specializes to this case as follows. First, one has to define a set of basis 
fields - naturally, for the static problem this is more easily done in terms of the potential rather 
than the field. At zero frequency, the Bloch wavenumber is also zero; hence the Bloch conditions 
for the field are periodic, but the potential may have an offset corresponding to the line integral of 
the field across the cell. Thus the first basis function ip\ is defined by 

V-eV^i=0; ^1 = ±|) = ±|; fo(y + a) - i/n(y) = 

The potential difference across the cell corresponds to a unit uniform applied field in the — re- 
direction. The second basis function is completely analogous and corresponds to a field in the 
— y-direction. Each of the basis functions can be found by expanding it into cylindrical harmonics 

Eoo 
c n g n (r, 4>); g n = (r n + p n r~ n ) sin(n0) 
n=l 

where the polar angle (f> is measured from the symmetry line of the potential (i.e. from the y axis for 
ijji and from the x axis for ^2); indexes 1 and 2 for the basis function and its respective expansion 
coefficients are dropped for simplicity of notation; the potential is gauged to zero at the center of 
the cell; p n is the polarizability of the inclusion for the nth harmonic. For a cylindrical particle, p n 
can immediately be found from the boundary conditions on its surface: 

. 1 ~ e cyl 2n 
Pcyl,n - 1 + ecy /cyl 

The potential at the x-boundaries is 



The coefficients c n can be found numerically by truncating the infinite series at re = n max and 
applying the Galerkin method. This leads to a system of equations 

Gc = f 
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Here c is the Euclidean vector of the expansion coefficients, 

G(n,m)= I g n g m dl; f n = \ [ g n dl 

Jx=a/2 Z Jx=a/2 

The system of equations may be solved numerically. 

Once the expansion coefficients and hence the basis functions have been found, the next steps 
are to compute the circulations and fluxes, and after that the WNBK interpolants, of the basis 
functions. For ipi, the circulation of the respective field ei = — V^i over each "horizontal" edge 
y = ±| is, by construction, equal to a and is zero over the two "vertical" edges. The fluxes of 
di = eei are zero over the horizontal edges; for the vertical ones, 



[[di]] = / d lx dy = - ed^idy 

Jy=±a/2 Jy=±a/2 



where e is that of the host if cell boundaries do not cut through the inclusions; the a>derivative can 
be computed from the r- and (^-derivatives of the cylindrical harmonics. 

Since the fluxes and circulations of the fields over the pairs of opposite edges are in this case 
equal, the WNBK interpolant is simply constant and the effective parameters are obtained imme- 
diately, without the intermediate stage of their pointwise values. Formally, the system of equations 
for e c g is 

'1 0\ fe xx e xy \ /[[di]] 
,0 l) \e yx eyy) { [[d 2 ]l 

the identity matrix is written out explicitly to clarify its origin: its first column represents the 
WNBK curl-conforming interpolant of ei; the second column - that of e 2 . The rightmost matrix 
contains the div-conforming WNBK interpolants of di and d2. As can be expected, the system of 
equations is in this case trivial. Since the fluxes [[di]] and [[d2]] are equal, the permittivity is, as 
expected, a scalar quantity equal to these fluxes. 

For numerical illustration, let us consider a cylindrical inclusion with a varying radius. The 
plots of e c ff vs radius of the cylinder (Fig. 5) for two different values of the permittivity (e cy i = 5 
and e cy i = 10) show an excellent agreement between the new method and M-G, in the range where 
such an agreement is to be expected - that is, for small radii of the inclusion. 



5.4 Example: Bloch Bands and Wave Refraction 

Definitive tests of the effective parameters come from Bloch band diagrams and, even more impor- 
tantly, from wave reflection and refraction at material interfaces. Here we consider wave propagation 
through a photonic crystal (PhC) slab - an array of cylindrical rods with no defects. For consis- 
tency with previous work, the geometric and physical parameters of the crystal are chosen to be 
the same as in the tests reported previously [32] and are taken from [9] . Namely, the radius of the 
rod is r cy i = 0.33a and its dielectric permittivity is e cy i = 9.61. The p-mode (ff-mode, with the 
one-component magnetic field along the rods) is considered because this is a more interesting case 
for homogenization. 

The numerical simulation of the Bloch bands and of wave propagation through the PhC slab 
was performed with FLAME [32, 33], a generalized finite difference calculus that is, arguably, the 
most accurate method for this type of problem because it incorporates local analytical solutions 
into the difference scheme. 
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0.1 0.2 0.3 

cylinder radius 

Figure 5: Effective e for a 2D-periodic array of cylinders. Lines: the proposed procedure (with 9 
cylindrical harmonics). Markers: the Maxwell-Garnett formula. 

Fig. 6 shows that the Bloch bands obtained with the effective parameters are in an excellent 
agreement with the accurate numerical simulation, except for a few data points at the band edge, 
where effective parameters cannot be expected to remain valid. 

Fig. 7 displays the real parts of e e g and /x e ff as a function of the Bloch wavenumber, in the TX 
direction. Among other features, a region of double-negative parameters (e c g < 0, /i e g < 0) can be 
clearly identified for the normalized Bloch wavenumber approximately between 0.33 and 0.42. This 
agrees very well with the band diagrams above and with the results reported previously [9, 32]. 

Numerical simulation (using FLAME) of EM waves in the actual PhC is compared with the 
analytical solution for a homogeneous slab of the same thickness as the PhC and with effective 
material parameters e c g and /i e g calculated as the methodology of this paper prescribes. In the 
numerical simulations, obviously, a finite array of cylindrical rods has to be used - in this case, 
24 x 5 to limit the computational cost. The results reported below are at the normalized frequency of 
Cj = Lua/(27rc) = (A/a)" 1 = 0.24959, coinciding with one of the data points in previous simulations 
[32]. A fragment of the contour plot of Ke(H) in Fig. 8 helps to visualize the wave for the angle of 
incidence tt/3. 

The real and imaginary parts of the numerical and analytical magnetic fields along the line 
perpendicular to the slab are plotted in Figs. 9 and 10, respectively, for normal incidence and in 
Figs. 11 and 12 for the angle of incidence tt/6. The analytical and numerical results are seen to 
be in a good agreement; some discrepancies can be explained by (i) numerical artifacts due to the 
finite length of the PhC in the tangential direction; (ii) the finite thickness of the slab that makes it 
different from the bulk; (iii) the approximate nature of the effective parameters, especially in this 
case where the cell size is ~l/4 of the wavelength and the filling factor is high. 
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Figure 6: The FX Bloch bands obtained with the effective parameters (markers) vs. accurate 
numerical simulation (solid line). e cy i = 9.61, r cy \ = 0.33a. 
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Figure 7: Effective parameters (markers - e' eS , solid line - n' eS ) vs the Bloch wavenumber in the 
YX direction. e cy i = 9.61, r cy i = 0.33a. 
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Contour plot of Re(H). 9i„ c — 7r/6, r ay i ~ 0.33, e cy i — 9.61 




Figure 8: Contour plot of Ke(H). Angle of incidence 7r/6, e cy i = 9.61, r cy i = 0.33a. 




Figure 9: Re(H) along the coordinate perpendicular to the slab. Normal incidence. Other param- 
eters same as above. As expected, the "microscopic" field (markers) exhibits some scatter within 
the slab, and nevertheless agrees fairly well with the effective medium field (solid line). The reasons 
for soem discrepancy are noted in the text. 
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Im(H) in the direction perpendicular to the slab. Normal incidence, r cy ; = 0.33, e cy i = 9.61 
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Figure 10: Im(H) along the coordinate perpendicular to the slab. Normal incidence. e cy \ = 9. 
r cy i/a = 0.33. 
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Re(H) in the direction perpendicular to the slab. 0; nt . = tt/6, r cy i = 0.33, e cy i = 9.61 
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Figure 11: Re(i/) along the coordinate perpendicular to the slab. The angle of incidence ir 
Other parameters same as above. 
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Figure 12: \m(H) along the coordinate perpendicular to the slab. The angle of incidence 7r/6. 
Other parameters same as above. 

5.5 Conclusion 

A new methodology is put forward for the evaluation of the effective parameters of electromagnetic 
and optical metamaterials. The main underlying principle is that the coarse-grained E and H fields 
have to be curl-conforming (that is, they have to possess a valid curl as a regular function, which 
in particular implies tangential continuity across material interfaces), while the B and D fields 
have to be div-conforming, with their normal components continuous across interface boundaries. 
While some flexibility in the choice of these coarse-grained fields exists, an excellent framework 
for their construction is provided by Whitney forms and the WNBK complex (Section 3.2). This 
construction ensures not only the proper continuity conditions for the respective fields, but also the 
compatibility of the respective interpolants, so that e.g. the curl of E lies in the same approximation 
space as B. As a result, remarkably, Maxwell's equations for the coarse-grained fields are satisfied 
exactly. 

Further, the electromagnetic field is approximated with a linear combination of basis functions, 
the most general choice of which is Bloch waves, although in special cases other options may 
be available. Effective parameters are then devised to provide the most accurate linear relation 
between the WNBK interpolants of the basis functions. In the limiting case of a vanishingly small 
cell size, this procedure yields the exact result; for a finite cell size - as must be the case for all 
metamaterials of interest [16, 34] - the effective parameters are an approximation. Ways to improve 
the accuracy are outlined in Section 4.3. 

Proponents of the differential-geometric treatment of electromagnetic theory have long argued 
that the H, E and B, D fields are actually different physical entities, the first pair best characterized 
via circulations (mathematically, as 1-forms) and the second one via fluxes (2-forms) [30, 12, 7]. 
The approach advocated and verified in this paper buttresses this viewpoint. 

There are several important issues to be addressed in future work. First, the theory of this 
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paper needs to be applied to 3D structures of practical interest. Second, physical considerations 
indicate that the theory should lead to causal and passive effective media (see Section 4.4), but 
a rigorous mathematical analysis or, in the absence of that, accumulated numerical evidence are 
highly desirable. The significance of passivity and causality is evident and has been emphasized by 
Simovski & Tretyakov [25, 26]. 

Although the object of interest in this paper is artificial metamaterials, it is hoped that the 
new theory will also help to understand more deeply the nature of the fields in natural materials 
because rigorous definitions of such fields, especially of the H field, are nontrivial. The ideas and 
methodology of the paper are general and should find applications beyond electromagnetism, for 
example in acoustics and elasticity. 
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